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04 . Abstract. 

We present strongly stable semi-discrete finite difference approximations to the 
quarter space problem {x > 0, t > 0) for the first order in time, second order in space 
wave equation with a shift term. We consider space-like (pure outflow) and time-like 
' boundaries, with either second or fourth order accuracy. These discrete boundary 

conditions suggest a general prescription for boundary conditions in finite difference 
I codes approximating first order in time, second order in space hyperbolic problems, 

' such as those that appear in numerical relativity. As an example we construct boundary 

conditions for the Nagy-Ortiz-Reula formulation of the Einstein equations coupled to 
a scalar field in spherical symmetry. 
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! 1. Introduction 

One of the major obstacles to obtaining long-term stability in numerical simulations of 
^ ■ strongly gravitating systems, such as a binary black hole space-time, is the proper 
■ treatment of boundaries. In general, when hyperbolic formulations of Einstein's 
equations based on space-hke hypersurfaces are used, one can distinguish between two 
types of boundaries: inner and outer boundaries. Whereas the inner boundary is a 
space-like hypersurface introduced to excise the singularity from the computational 
domain, the outer boundary is an artificial time-like surface introduced because of 
limited computational resources. 

If maximally dissipative boundary conditions are used with symmetrizable 
hyperbolic fully first order formulations of the Einstein equations El El El Ej it is 
known how to construct stable schemes to high order of accuracy. (Another desirable 
property of boundary conditions in general relativity is that they are compatible with the 
constraints, but we do not consider this here.) With first order in time and second order 
in space formulations such as [HI Ej , on the other hand, although the continuum problem 
is reasonably well understood |H1 El IHH EI] , niuch less is known about discretisations. 
The issue of constructing stable finite difference approximations of boundary conditions 
for first order in time and second order in space formulations of Einstein's equations has 
not yet been investigated and it is the focus of this work. As a first step we consider 
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the shifted wave equation in one space dimension and look at both second and fourth- 
order accuracy. A one-dimensional toy problem already captures many of the technical 
difficulties that arise in the higher dimensional case. For recent progress regarding fully 
first order or fully second order formulations, see OIISI- 

In Section |21 we review the continuum initial-boundary value problem for the 
shifted wave equation, both to fix the notation and to state the energy estimates that 
establish well-posedness of the continuum problem. Section El states our main results: a 
prescription for strongly stable finite differencing schemes for the quarter space problem 
(x > 0, t > 0), with second and fourth order accuracy, and for the two cases where the 
boundary is time-like or space-like (outflow). 

The main part of the paper is the proof of the stability and accuracy of these 
prescriptions at the semi-discrete level. Instead of using an energy method, we use the 
Laplace transform method, as described in Chapter 12 of ^1]. In its flnal step our proof 
relies on plotting a function of a complex variable to show that it is bounded. In SectionlH 
we apply this method to the semi-discrete initial-boundary value problem for the shifted 
wave equation and prove strong stability, estimate ()84p . and convergence. Finally, we 
turn the semi-discrete scheme into a fully discrete one using the fourth-order Runge- 
Kutta method for integrating in time. In SectionElwe present numerical tests confirming 
the desired order of convergence of the resulting schemes for the wave equation. As an 
application, in Section IHl we consider the Nagy-Ortiz-Reula formulation of the Einstein 
equations in spherical symmetry, where the boundary conditions introduced in Sectional 
are generalized to constraint-preserving boundary conditions. 

For the second-order in space wave equation without a shift term, stable and 
accurate boundary conditions can be constructed using the summation by parts rule 
[T31 IT^ . The proof of stability and accuracy of these methods uses the existence of a 
discrete conserved energy which is the precise analogue of the continuum energy. In 
[Appendix A| we attempt to generalise these methods to the wave equation with a shift, 
but we fail. The reason is that in the presence of a shift, three separate summation by 
parts properties must be obeyed instead of one with vanishing shift, and this appears 
to overspecify the finite differencing scheme. This suggests to us that standard discrete 
energy methods are not suitable for the second-order in space wave equation with a 
shift, and by extension for other second-order in space hyperbolic systems. 

For completeness, and for comparison with the wave equation, we give the analogous 
results for the advection equation in [Appendix B In [Appendix C we briefly compare 
our results with those of 



2. The continuum initial-boundary value problem 

The one-dimensional shifted wave equation consists of a system of coupled linear partial 
differential equations (PDEs) of the form 

dt(j) = (3d,(j) + n + (1) 

dtU = (3d,U + dl<P + F", (2) 
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where F'^{x, t) and F^{x, t) are forcing terms and the parameter (3, the shift, is assumed 
to be constant. The homogeneous version of system (PI21) is obtained from the one- 
dimensional wave equation, d?(j) = after a Gahlean change of coordinates, t = i, 
X = X — pi and the introduction of a new variable 11 = dt4> — Pd^cp- (See [T7] and 
[Appendix C| for why this is better than using dtcj) in finite difference schemes.) At time 
t = we prescribe initial data 

0(x,O) =/<^(x), (3) 
n(x,0) = /n(x). (4) 

Introducing the vector valued function u{x,t) = {(f){x,t),Il{x,t),dx4>{x,t))^ , with 
periodic boundary conditions on an interval V the energy 



[ [<p' + u' + {d,<py]dx (5) 

Jv 



satisfies the estimate 

lk(-,t)f <^w(lk(-,0)||2 + ^*l|F(-,r)fdr), (6) 

where F{x, t) = (F'^(a;, t), F^{x, t), d^F^^x, t))^ and K{t) is a function which is bounded 
on any finite time interval and does not depend on the initial data. 

We now introduce a boundary at a; = and consider two "quarter space" problems 
(meaning x > and t > 0): one for the pure outflow case (/3 > 1) and one for the 
time-like case (|/5| < 1). No boundary condition is needed in the outflow case. In the 
time-like case we impose the Sommerfeld boundary condition 

w.(0,t) = (7(t), (7) 

where w± = 11 ±9x0 are the characteristic variables and g is a, freely specifiable function 
compatible with the initial data. To obtain an energy estimate we take a time derivative 
of the energy and use integration by parts to obtain the inequality 

^||w(-, t)f< -\pi[^ + 2n9,.0 + P{d,(t^)%=o + const {\\u{; t) f + ||F(-, t) f) . (8) 
Rewriting the boundary term as 

i[(l-/3)^!-(l+/3X].=o, (9) 

which is negative definite in the outflow case and bounded by |(1 — P)g'^ in the time-like 
case, and integrating, gives the estimate showing strong well-posedness 

\H-M? < m (ll/f + j\\\F{., r)f + 6\g{r)\')dr^ , (10) 

where f{x) = u{x, 0) and 5 = for the outflow case and 5 = 1 in the time-like case. 

Using the energy method we have proved the well-known fact that the initial- 
boundary value problem for the shifted wave equation is well-posed in the outflow case 
with no boundary condition, and in the time-like boundary case with a Sommerfeld 
boundary condition. In the remainder of this paper we investigate the stability of 
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finite difference discretisations of tliese two problems. We sliall restrict our discussions 
to discretising tfie equations in space but not in time (tlie metliod of lines). Tfiis 
transforms tlie partial differential equation into a large coupled system of ordinary 
differential equations (the semi-discrete problem) which can be solved by a standard 
ODE integrator, such as fourth order Runge-Kutta. 



3. The semi-discrete initial-boundary value problem: summary of results 

In this Section we summarize our results regarding strong stability and convergence 
of schemes approximating the quarter space problem for the shifted wave equation. 
We introduce the grid xj = jh, with j = 0, 1,2, . . ., and the grid functions (f)j{t) and 
Ilj{t) approximating the continuum variables. In the interior, j > p/2 (with p = 2,4 
depending on the accuracy of the scheme), we use the standard centered minimal width 
discretization 

^0^- = + n, + F,^ (11) 
An, = /?D«n, + z}(2)0, + i^n (12) 

where D^^'^ and D*-^-* approximate the first and second derivatives, respectively. In the 
second order accurate case, these operators are 

D^^^ = Do, = D+D_, (13) 

where D^Uj = {uj^i — Uj)/h, D_Uj = {uj — Uj_i)/h and Dq = [Dj^ + DJ)/2. In the 
fourth order accurate case we use 

= D,{l- ^D^D^j , (14) 
= D+D^ ( 1 - —D+d\ . (15) 



12 

We know that the Cauchy problem (no boundaries) for ()11H12|) is stable for any 
value of the /? in both second and fourth-order accuracy |.17^. 

For the quarter space problem, it is convenient to introduce ghost points, that is, we 
assume that the interior equations hold for all j > and provide numerical prescriptions 
for the grid functions at the grid points j = —p/2, . . . , — 1. Stable discrete boundary 
conditions are given in the following subsections. The proofs follow in Section |3J 

For completeness, and for comparison with the second-order wave equation, we give 
second and fourth-order accurate boundary prescriptions for the advection equation in 
Appendix B 



3.1. Second order accuracy 



3.1.1. Outflow boundary We start with the outflow case (3 > 1. The continuum 
problem does not require any boundary condition, but at the discrete level a special 
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prescription is needed. Third order extrapolation for and second order extrapolation 
for Ilj, namely 

= 0, (16) 

h'^DlU^i = 0, (17) 

give strong stability and second order convergence. Interestingly, the minimum order 
of extrapolation required for second order convergence is not the same for the two grid 
functions (pj and 11^. The reason for this becomes clear in the convergence analysis 
of subsection 14.1.11 This is to be contrasted with the result for fully first order 
symmetrizable hyperbolic systems, in which second order extrapolation (or, equivalently, 
first order one-sided differencing) for the outgoing characteristic variables yields second 
order convergence. 

3.1.2. Time-like boundary When the shift satisfies < 1, one of the two characteristic 
variables is entering the domain through the boundary. We seek discrete boundary 
conditions approximating 

n(O,t)-a,0(O,t) = ^7(t), (18) 

which lead to strong stability and preserve the internal accuracy. This is achieved by 
populating the ghost point j = — 1 using 

Ho - 1^000 = 9, (19) 

h^Dln^i = 0, (20) 

or, explicitly, 

0-1 =0i + 2/i((7-no), (21) 

n_i = 2Uo - Hi. (22) 

3.2. Fourth order accuracy 

3.2.1. Outflow boundary In the outflow case, the extrapolation conditions 

h^Dl(f)_i = 0, (23) 
h'Dl<f)^2 = 0, (24) 
h^DlU_i = 0, (25) 
/iXn-2 = 0, (26) 
lead to strong stability and fourth order convergence. 

3.2.2. Time-like boundary The prescriptions 

Ho - 1^(1)00 = g, (27) 

/i^Z}^0„2 = 0, (28) 

/i^D^n_i = 0, (29) 

h^DlU^2 = 0, (30) 
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where D^^'' is defined in (jl4j) . give strong stability and fourth order convergence. 
(Increasing the extrapolation order by one in both (p and 11 also gives fourth order 
convergence.) Explicitly solving ()27|1 and ()28j) for </)_i and 0_2 gives 

0_i = A{g - Uo)h - y0o + 601 - 202 + ^h, (31) 
0_2 = 20{g - Uo)h - y0o + 4001 - 1502 + ^03- (32) 

4. Proofs of strong stability and convergence 

In this Section we use the Laplace transform method for difference approximations 
as described in Chapter 12 of ^1] to prove strong stability for second and fourth- 
order accurate discretisations of the initial-boundary value problem for the shifted wave 
equation. In order to apply the theorems of that reference, which assume that hyperbolic 
problems are written in fully first order form (see equation (12.1.11) of Jl]), we need 
to perform a discrete reduction to first order [18j. 

4-1- Second order accuracy 

4-. 1.1. Outflow boundary We consider the semi-discrete quarter space problem for the 



outflow case (/? > 1) 

^0,- = PDo<P, + Uj + Ff, (33) 

^n, = pDoU^ + D+D_4>j + Ff, (34) 

0.(0) = //, (35) 

n,(0) = /f , (36) 

= g^, (37) 

/i^iL>^^n_i = g^, (38) 

l|n||^+P+0||^<oo, (39) 



where j = 0,1,2,..., \\u\\1 = Yl'jLo I'^iP^^ qi and q2 are non-negative integers. In 
practice one would choose = = 0, but in the analysis that follows we will need 
the inhomogeneous case. A first order reduction of the problem obtained by introducing 



the grid function Xj = D^(j)j gives 

^n, = pDoU, + D_X, + Ff, (40) 

^X,=/5DoX, + D+n, + if , (41) 

n,(0) = ff, (42) 

^.(0) = ff, (43) 

h'^W^U^i = g^, (44) 

/^<?2^92j^_^ = g^, (45) 
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where j = 0,1,2,.. ., Ff = D+Pf, ff = D+ff, and = g'^/h. The auxihary 
constraint Cj = Xj — D^(j)j satisfies the homogeneous system of ODEs 

^Cj = f3DoC,, J = 0,1,2,..., (46) 

cm = 0, (47) 
h'i^DfC.i = 0, (48) 

and therefore vanishes identically. 

We want to show that the semi-discrete initial-boundary value problem ()40M5|) is 
strongly stable and second order convergent if gi,g2 > 2. Since the evolution equation 
for (f)j only involves lower order terms {DQ(j)j can be expressed as a combination of Xj 
and Xj-i), it can be ignored in the analysis that follows. 

The proof of stability is divided into three steps. We first estimate the solution 
of the problem with F = and / = in terms of the boundary data by checking 
that the Kreiss condition, inequality (jM}, is satisfied. We then estimate the solution 
of an auxiliary problem with modified homogeneous boundary conditions. Finally, we 
combine the two estimates to obtain the strong stability estimate, inequality (jHSl), for 
the original problem. In its final step the proof relies on plotting a function of a complex 
variable to show that it is bounded. 

Step 1. By estimating the solution of the F = 0, f = problem near the boundary 
using the Kreiss condition ()64|) we obtain the estimate ()49|). Let F = and / = in 
Eqs. fl40ll45|l . We want to show that 

l|n(t)||^ + ||X(t)||^ < const fmr)\' + \g^{T)\') dr. (49) 
Observing that 

^(||n(t)||^ + lix(t)ii^) = -/5(non„i + XoX_o - 2noX_i (50) 



< const (in^p + iX^n, 
i=-i 

it is clear that if we can estimate the solution near the boundary (i.e., at j = —1,0) in 
terms of the boundary data, we recover the estimate ()49|). We show that this is possible 
by explicitly solving the Laplace transformed problem 

sfij = /?(n,+i - n,-i)/2 + X, - x,_i, (51) 

sXj = /3(X,+i - X,-_i)/2 + fl,-+i - n,, (52) 

/i^iD^^n.i = (53) 

h'^WfX., = (54) 
||flf,+ ||X||^^<oo, (55) 

where u{s) = f^°° e"'^*u(t)dt and s = sh. 
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Eqs. (15 ip and (|52j) form a system of difference equations. The characteristic 
equation associated with it, obtained by looking for solutions of the form flj = k^flo, 
Xj = yXo, is 



(56) 



a polynomial of degree 4 in k. For Re(s) > there are no solutions with \k\ = 1. If 
k = e'^ with ^ G M was a solution we would have Re(s) = 0, which is a contradiction. 
Observing that the roots are continuous functions of s and for large values of \s\ we 
have s ~ 13 /{2k) ± we conclude that for Re(s) > there are two and only two 

solutions, ki and /c2, inside the unit circle. For s = the four roots are 



. i^i?!±|v^. ,,,,.1. (57) 

Since, for small ^3 4 = 1 + (/? ± l)^^s + 0(|sp), the roots ki and k2 are those which 
are inside the unit circle for Re (5) > 0. 

The general solution of the difference equation ()5HI52|) satisfying ||n||^ + < 00 
can be written as 

lij = aisik{ + (r2S2ki, (58) 
X. = a,{ki - l)ki + a2{k2 - l)ki, (59) 

where §1,2 = s — f (^^1,2 — ^1^2) ■ Note that since we will be constructing the explicit 
solution, the case ki = k2 does not require special treatment. Inserting (|58H59|) into the 
boundary conditions ()53II54|) gives rise to a 2 x 2 system in o"i and o"2 

a^hih - ly'K^ + (y2Hk2 - ifK^ = g^, (60) 

ai{ki - ly^+^k^' + a2{k2 - lY'+^k2^ = g"". (61) 

If the coefficient matrix C{s) multiplying (ai, (J2Y' is non-singular for Re(s) > 0, we can 
solve for ai and (J2 and substitute into and (jKHjl . and obtain an explicit solution of 
the Laplace transformed problem, which we write in the form 

n. = E 4^'' (62) 

fc=n,x 

k=Il,X 

To verify the Kreiss condition, namely that 

x:(in.r+i^.r)<^(i^"r+irr), (64) 

i=-i 

we numerically compute the coefficients c^^, c^, for j = —1,0 and A; = X, 11, and plot 
the quantity 



N 



E (i-y.P + i4i^) 



(65) 
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We restrict our attention to the qi = q2 = 2 case (same order of extrapolation for 11^ and 
Xj). Since for |s| > Cq we have that ifr^f + < Ko{\g^\'^ + \g^\'^), see Lemma 12.2.2 
of [HI, we only need to consider the compact set S* = {s G C : |s| < Cq, Re(s) > 0}. 
From the chain of inequalities 



||fl||2 + ||X||^ = ^ [\\hf3DoIl + hD_X\\l + \\hpDoX + hD+n\\l 
- ^ {^^\\^^''^\\h+ \\hD.X\\l + p^\\hDoX\\l+ \\hD+U\\l 



2 

< 



(3^ [ \\U\\l + -|n_i|2/i ) + 4||X||2 + 2|X_i|2/i + (3^ ( \\X\\l + ) + 4||n| 



<^,{^ + P'){ml + \\X\\l + \9''\'h+\g''\'h 



where in the last inequality we used the boundary conditions (1531154^ . we see that for 
P = 2 we can take Cq = 12. In figure [T] we plot N{s) where s E S. 




Figure 1. On the left we plot N{s) for s E S O {Im(s) > 0}. Since this function is 
symmetric across the real axis, we only need to display the region with non-negative 
imaginary part. The Kreiss condition is satisfied for the initial-boundary value problem 
(|33I39|) with /3 — 2 and qi = q2 = 2. The spike that appears at a point near Re(s) = 
is bounded. This is illustrated on the right, where we have increased the resolution of 
the plot in a neighbourhood of this point. The function N(s) is continuous but not 
differentiable on Re(s) = due to branch cuts in its definition. 



Using Parseval's relation and the fact that the solution at time ti does not depend 
on the boundary data at time t2 > ti, the Kreiss condition implies the estimate in 
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physical space 

r (|n,f + |X,n dr<K [\\gY + \g^\') dr. (66) 
^0 ^■^„i Jo 

Combining this with the integral of inequality (jKn|) gives the desired estimate fl49|l . 

Step 2: We estimate the solution of a problem with modified homogeneous boundary 
conditions in terms of the initial data and forcing term, inequalities ()74|) and (fTBj) . So 
far we have shown that, for vanishing initial data and in the absence of a forcing term, 
the solution can be estimated in terms of the boundary data. We now consider the 
auxiliary problem 

^n, = pDoU, + D_X, + Fj", (67) 

^X, = f3DoX, + D+U, + Ff, (68) 
n,(0) = /f, X,(0) = /f, (69) 

n_i = l(llo- Ixo] , (70) 



X-i = |xo, (71) 

||n||2 + ||X||2 <oo, (72) 

where the boundary conditions were chosen so that a direct estimate in physical space 
can be obtained. The estimate 

^(linil^ + \\x\\l)< -2(|nop + iXoH + ||n||^ + ||x||^ + iif^h^ + (73) 

implies 

l|n(t)||^ + ||X(t)||^ < const (^WfWl + £ llF(r)ll^dr) (74) 



and 

ft 



(75) 



^ (|n,f + |X,p) dr < const (^\\f^\\l + ll/'^ll^ + ^ dl^^ll' + WF^'WDdr 

for j = —1,0. Since our interior scheme uses the same number of grid points in each 
side, one can show (Lemma 12.2.10 of 14 ) that for every fixed j 

POO / POO \ 

I (|n,r + |X,f)dt< const (^11/11^ + \\F{t)\\ldtj , (76) 
where / = (/", f^f and F = (F", F^f. 



Step 3: Using the estimates of Steps 1 and 2 we derive the estimate ()83j) showing strong 
stability. If we denote by 11^ and the solution of the auxiliary problem of Step 
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2, and by 11^ and Xj the solution of the original problem, we see that the differences 
f[j = Uj - n| and Xj = Xj - X^ satisfy 

jn, = pDon,+D.X^, (77) 

^^X,=f3DoXj + D+nj, (78) 

n,(0)=0, (79) 

X,(0)=0, (80) 

h'W_^fl^i = g^- /i^^D^^n^i, (81) 

h'^WfX^i = g^ - h'^WfX^. (82) 
Using ()76|) and (f?^ we obtain the estimate 

l|n(t)ll^ + ||x(t)||^, < 2(||n(t)||^, + ||x(t)||^, + ||n^(t)||^, + ||x^(t)||^,) < (83) 

const (||/||^ + ^\||F(r)||^ + |^(r)r) dr^ 

where g = {g^,g^)'^. This inequality proves strong stability. Reintroducing the 
evolution equation for (pj, we have an estimate with respect to the D+-norm 

+ l|n(t)||^, + WD^mWl < const(||/^||^ + Wf^Wl + \\D^f% + (84) 
{\\FHr)\\l + ||i^"(r)||^ + \\D^FHr)\\l + \g^{r)\' + |^<^(r)//^ndr) . 



Proof of convergence: Having shown strong stability it is straightforward to prove 
convergence. We only need an estimate for the error. Defining the errors e^(t) = 
(pjit) — (j){xj,t), e^{t) = Ilj{t) — Il{xj,t) we see that they satisfy the initial-boundary 
value problem 

= (3D,et + ef + Ff, (85) 

^ef = PDoef + D+D^e'^ + Ff, (86) 

e,^(0) = 0, (87) 
6^(0) = 0, (88) 
/^<72+i^a2+ig^^ = -/,'?2+i0fe+i)(o, t) + 0(/i'?^+2), (89) 

h'iW^e\ = -/i'?in(^i)(0,t) + 0(/i^^+^), (90) 

where Ff = 0(/i^) and Fj^ = 0(/i^). We perform a discrete reduction by introducing 
the quantity ef = -D+e^ which satisfies 

^ef = /3Doef + D^e"} + F^ , (91) 

ef (0) = 0, (92) 
h'^^Dfe^ = -/i'?2 0('?2+^)(O,t) + 0{h''^+^). (93) 
Note that since F/ = /3{Do(l)ixj,t) - (f)^{xj,t)), we have Ff = D+Ff = Oih"^). 
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The strong stability estimate (|S^ applied to (|SH|). (|SS|). (j9UH93p . guarantees that 
the scheme is second-order convergent, i.e. 

(l|e1l^ + ||e"||^+P+e^||D'/'<0(/^^), (94) 

provided that qi,q2 > 2. In this analysis we have implicitly assumed that the initial 
data and forcing terms are exact. However, this assumption can be easily replaced by 
the requirement that the initial data errors for and ef are of order h"^ (note that this 
means that Dj^e- = 0(/;,^)) and that the errors in the forcing terms Fj^ and Fj^ are of 
order h"^. An immediate consequence of (jMj) is that we have convergence with respect 
to the discrete L2-norm, (||0||^ + ||n||^)^/^. 

We have also studied the case qi = q2 + 1 = 3 (same order of extrapolation for (pj 
and Uj). The Kreiss condition can be verified directly and the rest of the stability proof 
applies. 



4.1.2. Time-like boundary In this Subsection we prove stability for the shifted wave 
equation problem with a discretization of the Sommerfeld boundary condition, equation 
()18|). The case with zero shift was discussed in Appendix A of using the discrete 
energy method. Here we focus on the 7^ \(3\ < 1 case. 

We consider the same semi-discrete evolution system (jHHHHfjj) and with 
boundary conditions 

Ho - Do0o = g"", (95) 
hWlU_, = g^. (96) 

In applications one would set g'^ equal to the boundary data g{t) of the continuum 
problem and g^ = 0. The discrete reduction of is 

Ho-i(Xo + X_i)=^7'', (97) 

which implies the boundary condition C_i = —Cq for the auxiliary constraint. 

The proof of strong stability proceeds as in the outflow case. We only need to show 
that the Kreiss condition is satisfied. Inserting Ilj{t) = e^^/c-'Ho and Xj{t) = e^^k^X^ 
into the scheme, and looking for non trivial solutions gives the characteristic equation 
fj56|) . For small |s| the roots are 



t.- '-'''-y^ ,Om, (98) 

fc2 = l + -^5 + 0(|sT), (99) 
,,.l^£±l^,Om, (100) 

*4 = l + ^3 + 0(|st). (101) 

For Re(s) > the roots ki and k2 have magnitude smaller than 1 and the remaining 
two have magnitude greater than 1. The requirement ||n||^ + \\X\\1 < 00 implies that 
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the general solution will have the form (j58|) and (j59|) . where the parameters cxi and (T2 
can be determined by imposing the boundary conditions 

cri - i(l + mi - K')) - y(1 + (1 + ^) = 9^ (102) 

aMh - irK' + (^2(A:2 - irk,' - ~si{h - l^K') = (103) 

The solution, if it exists, has the form (|62II63|) . Again, we verify the Kreiss condition 
by plotting the quantity for 5 G 5* and verifying that is is bounded. 

Having established that the estimate ()49|1 holds, we use the same auxiliary problem 
used in the outflow case. Step 2, giving the estimates (f7i|) and (17^]) . Hence, we have 
strong stability. Convergence follows by observing that the error equation associated 
with (jHSD is 

e^-Doet = ^4>"\0,t) + 0{h'). (104) 


4-2. Fourth order accuracy 

4-2.1. Outflow boundary The fourth order accurate standard discretization of the 
shifted wave equation is given by 

^0,- = + n, + Ff, (105) 

^H, = /?D(i)H, + D^^Uj + Ff, (106) 
at 

where j = 0, 1, 2, . . . and 

D^'\^ = Do(^- ^D+D^ Uj = + 8m,+i - 8m,_i + Mj_2)/(12/i), (107) 

D^'^^Uj = D+D^ (^1 - Y^D+D_^ Uj (108) 

= (—Uj+2 + I6uj+i — 30Mj + 16Mj_i — 'Uj-2)/ (12/i^). 
We consider the outflow case first [fi > 1), with the boundary conditions 

h'Dl(l)^, = gt, h'Dl<l).2 = 9t, (109) 
h'DlU^,=gY, h'DlU^2 = 9f. (110) 
As in the second order accurate case, we need to perform a discrete reduction to 
first order. For this purpose it is convenient to introduce a grid function Xj, which 
is a suitable linear combination of D^(j)j and D_(f)j. The choice of the particular 
combination is determined by the following two observations. First, the operator D^'^^ 
can be decomposed as 

D^'^^ = D+D^, (111) 

where D± = (1 =1= ahD^)D± = (1 — a)D± + aD^ and a is a root of the quadratic 
equation — a — 1/12 = 0. Second, in the absence of boundaries, the following discrete 
energy, 

imi+WD^m + ^WD+D^m (112) 
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which is equivalent to ||n||| + is conserved and it can be written asj 

l|n||^ + ||x||2, (113) 

where Xj = D^(j)j. 

These two facts suggest introducing the discrete variable Xj, leading to the 
(interior) discrete reduction 

^n, =pD^'% + D^Xj + Fj', (114) 

^X, = pD^'^Xj + D+U, + , (115) 

n,(0) = ff, (116) 

^.(0) = ff, (117) 

where Ff = D+Ff ff = D+ff. 

We need to translate the extrapolation boundary conditions for (p and 11 in terms 
of boundary conditions for the new variables, 11 and X. To do this, we go back to 
the original system ()l()5lll()(ij) and eliminate the variables 0_2, n_i, and n_2 using 
the boundary conditions ()109II110|) . Defining Xj = D+(f)j for j > as before, where in 
Xq the grid function 0_i is eliminated using ()l()9j) but with gf = 0, we can eliminate 
each occurrence of (pj in terms of Xj. The evolution equations for the Xj variables can 
be computed by taking appropriate combinations of the evolution equations of the (pj 
variables. This leads to a semi-discrete problem for Ilj and Xj, j > 0. For j > 2 the 
equations have the form ()114II115|) . with the exception that 

Ff = D+Fi + aPgt/{12h'). (118) 

However, for j = 0, 1 they are more complicated. To analyze the stability of the system 
we reintroduce ghost points. By setting the evolution equations near the boundary to 
be formally equal those at the interior, for j = 0, 1, we obtain the prescriptions 

n_i = 4no - eni + 402 - + + ^, (119) 

n_2 = lOHo - 20ni + ISHa - 4n3 + 4(?n + + ^^Lz^^^f , (120) 

X_i = 4X0 - 6X1 + 4X2 - X3 + (1 - 5a)g^ + ag^, (121) 
_ -137 + 132a -289 + 288a -241 + 252a ^ -43 + 48a ^^ 

^ 12a - 13 ° 12a - 13 ^ 12a - 13 ^ 12a - 13 ^ ^ ' 

12a - 11 ~ 1 ~ a - 1 ^ a - 1 ^ 
-X4 + — — X5 + 12—— —no - 48—— —Hi 



12a -13 12a -13 /5(12a - 13) /3(12a-13)" 

+72 ^ na - 48 ^ + 12 ^ ^ U^ 

/3(12a-13) /?(12a-13) ^ /?(12a - 13) 

2(6 - 6a + 53a/?2 _ 55/^2) ^ a-2^ a-1 „ 

+ /32(12a - 13) + 12a - 13^^ - 12^(i2a - 13)^^ ' 

i To show that \\D+<j)j - ahD+D^(f)j\\l = \\D+(j}\\l + /12\\D+D_(t)\\l one can use the identity 
hD+D^ = D+-D^ and _ ^ - 1/12 = 0. 
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where = gf/h. The forcing terms near the boundary are = D+Ff and 
= b^F^ (assuming Ft^ to be defined via h^DlFti = 0). 

We now need to show that the semi-discrete initial-boundary value problem ()114l - 
I122j) is strongly stable and fourth order convergent. Neglecting forcing terms, we have 
the estimate 

^(||n||2 + \\x\\l) = ^(x_iXi + n_ini + x^^Xo + u^^Ho - 8UoU_, - sx^M 

1 

+2{a - l)X_ino - 2an_iXo < const ^ (jn^f + \Xj\^). 

i=-2 

As in the second order accurate case, we explicitly solve the Laplace transformed problem 
for vanishing initial data and no forcing term and write the solution as 

n, = c%gf + c%g^ + c^^^^f + c%^g^, (123) 
X, = cfn.^r + cf^M + cfxM + cfx,92- (124) 
We numerically verify the Kreiss condition by plotting the quantity 

N= E + ■ (125) 

\ k = ni,n'2,x[,X2 ) 

inside the semi-disk S with Cq = 30. 

The modified homogeneous boundary conditions for the auxiliary problem are 

X_i = (^-48Xo - 6Xi + ^^^y ^^ Ho^ , (126) 
X_2 = - 4(65/5' + 144a2)Xo, (127) 

n_i = (48^no - en, - ^Xo) , (128) 

n_2 = - -r^((65a - + 12(a - l))no - —Hi, (129) 

p^a pa 

and they give the estimate 

^(linil^ + llxilD < - E(|n/ + 1^/) + linil^ + + llF^n^ + (130) 

i=o 

This implies inequalities fl74|) and (f7B|) and strong stability with respect to the norm 
follows from the fact that ||n||| + ||X||^is equivalent to ||n||^ + ||_D+0||^. Convergence 
is a consequence of estimates for the error. However, due to the modified forcing term 
(jllSp we are only able to show that ||n||^ + ||X||^ < 0(/i^). Given that the numerical 
tests of Section El indicate fourth order convergence, we believe that this estimate is not 
optimal. 
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4-2.2. Time-like boundary In the time-like case (0 ^ \(3\ < 1) the discrete boundary 
conditions are 

Ho - D^^Vo = gf, (131) 
h'Dl<f).2 =9t, (132) 
/^Xn-i =9f, (133) 
h'DlU_, = g^, (134) 

where D^^^ is given in (fT^ . We introduce Xj = D+cpj for j > 0, where 0_i in Xq is 
given by system ()131II132|) with = 0. Note that Xq contains also Ho. A discrete 
reduction to first order gives the ghost-point prescriptions 

n_i = 4no - 6ni + -U^ + gf- ^g^ + ^g^, (135) 

n_2 = lOHo - 20ni + 15n2 - + 4^" + g^ - 32^^^^-^g^ (136) 

8 16a - 17 X 
^3/3(10a-9)^2 ' 

X_, = -^^Xo + l^X,-^^X.-^^no-^^2^,f (137) 
10a - 9 10a - 9 10a - 9 10a - 9 3 10a - 9 ^ ^ 

8 3a - 4 y 

H 99 , 

9 10a -9^^ ' 

~ _ -383/3 - 1556 + 1448a + 358/3a ~ 2 528a - 552 + 866/3a - 1001/? ~ 

- /?(118a-127) ^° + 3 /5(118a-127) ^^^^^ 

_1690/3a+ 168a - 180 - 773/?^ 2 54a - 59 ^ 1 10a - 9 ^ 
~3 /?(118a - 127) ^ ^ 3 118a - 127 ^ ~ 3 118a - 127 ^ 

4 666/5a - 287 - 697/5 + 270a ^ _ ^^59 + 54/5a - 54a - 59/?^ 
^3 /3(118a- 127) °~ /?(118a - 127) ^ 

^^59 -I- 54/5a - 54a - 59/5^ 8 59 + 54/5a - 54a - 59/3^ 
^ /3(-127+ 118a) ^~3 /5(118a - 127) ^ 
_2 708 + 2376a/5 - 2596/? - 4051/3^ - 648a + 3750ap^ x 
~3 /52(118a - 127) 

-648a + 1728a/5 - 5865/3^ + 708 - 1888/3 + 5442a/32 ^ 
^ 18/32(118a- 127) ^2 

2 -162a + 216/3a- 236/3 + 177 „ 2 54a - 59 „ 
^3 /3(118a - 127) ~ 3118a - 127^^ ' 

where g^ = gf and g2 = gt/h. The forcing terms near the boundary are 

Fi^ = D^Fi + ^^ig^-12g^), 
= D^Ff, 

and F^ is obtained from the definition of Xq with the replacements (pj Fj^, 11^ Fj^. 

The rest of the proof proceeds as in the outflow case. That the scheme is convergent 
follows from estimates for the errors. As in the outflow case, we are only able to show 
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3.5th order convergence, although experiments suggest that the scheme is fourth order 
convergent. 

5. The fully discrete initial-boundary value problem: numerical tests 

The semi-discrete schemes we have considered here can be turned into fully discrete 
schemes using, for example, fourth order Runge-Kutta as a time integrator. In general, 
stability of the semi-discrete scheme does not guarantee that the fully discrete scheme 
is stable |2II1 E] • We therefore check stability of the fully discrete scheme by performing 
numerical convergence tests. 

The convergence tests are performed using the domain < x < L, where x = is 
the physical boundary, and by monitoring the errors in the interval < x < 1 at time 
t = 1. The value of L is chosen large enough so that the numerical solution in the interval 
< X < 1 for t < 1 is numerically unaffected by the boundary at x = L. We have used 
a Courant factor of 0.5 and resolutions ranging from h = 1/25 to 1/400. The code is 
tested against the exact solution 0(x, t) = f{—x + (1 — n(x, t) = /'(— x + (1 — P)t), 
where /(x) = sin(27rx). 

As shown in table ^ we find good second or fourth-order convergence in the norm 
(||n||^ + ||Z}+0||^)^/^ over all resolutions. We focused on the f3 = 2 and /? = — 1/5 cases, 
but the schemes are convergent for other values of the shift parameter in the ranges 
(3 > 1 and < 1. However, we noticed that in the time- like, fourth order accurate 
case, in order to avoid losing accuracy, the boundary data g had to be Taylor expanded 
at the intermediate time steps of the Runge-Kutta integrator. This issue is discussed in 

m- 

6. Application to the Nagy-Ortiz-Reula formulation of the Einstein 
equations in spherical symmetry 

The main interest of the one dimensional shifted wave equation IS cLS cL toy model 
for second order in space formulations of three-dimensional general relativistic initial- 
boundary value problems. For such problems initial and boundary data cannot be freely 
specified. The initial data has to satisfy the constraints and the boundary conditions 
should be such that no constraint violation is injected. Whereas maximally dissipative 
boundary conditions involve only first derivatives across the boundary, constraint- 
preserving boundary conditions require second derivatives across the boundary. 

We consider the Nagy-Ortiz-Reula (NOR) formulation 23^ of the Einstein 
equations, restricted to spherical symmetry and coupled to a massless minimally coupled 
scalar field. In spherical symmetry the line element takes the form 

ds'^ = -a^dt^ + grr{dr + f3dtf + geeidO"^ + sin^ Odcf). (139) 

and the reduction of the NOR system to spherical symmetry is straightforward. 
However, if the origin is included in the domain, the definition of the auxiliary variable 
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Table 1. The discrete boundary conditions described in Section 01 give the expected 
order of convergence. 

Errors and convergence rates 
T2 q 



(a) Second order accurate case, /3 = 2 



25 


7.3508410- 


-1 




50 


1.8395110- 


-1 


1.9986 


100 


4.60081 10- 


-2 


1.9994 


200 


1.1502110- 


-2 


2.0000 


400 


2.8755510- 


-3 


2.0000 


(b) 


Second order accurate 


case, P = 


25 


1.0604210- 


-1 




50 


2.59231 10- 


-2 


2.0323 


100 


6.4155910- 


-3 


2.0146 


200 


1.5967310" 


-3 


2.0065 


400 


3.9836610" 


-4 


2.0030 



-1/5 



(c) Fourth order accurate case, P = 2 



25 


9.7074710" 


-3 




50 


6.1033410" 


-4 


3.9914 


100 


3.8202410" 


-5 


3.9979 


200 


2.3880910" 


-6 


3.9997 


400 


1.49255 10" 


-7 


4.0000 


(d) 


Fourth order accurate 


case, P = 


25 


1.0195510" 


-3 




50 


6.7179010" 


-5 


3.9238 


100 


4.3012810" 


-6 


3.9652 


200 


2.7206410" 


-7 


3.9827 


400 


1.7102510" 


-8 


3.9917 



-1/5 



/ and the addition of its definition constraint to the right hand side of the evolution 
equations require consideration. We use the regularised dynamical variables gx = r^'^gee, 
Kt = r~'^Kgg, and / = —2(\ngT),r, and instead of adding a multiple of Q^r to the right- 
hand side of the evolution equation for Krr, where Q = f + 2{\n gT),r is the definition 
constraint for the auxiliary variable /, we add Q^r~Q / f- The regularity conditions are the 
origin are that grr, gr, Krr, Kt, 4>, H are even functions of r satisfying gT\r=Q = grr\r=o, 
KT\r=o = Krr\r=o and / is an odd function of r. 

With fixed densitised lapse, a = y/g^grQ, and fixed shift, the evolution equations 

are 

dtgrr = Pdrgrr + 2grrdrP — 2aKrr, (140) 
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(141) 

drgrdr.Q (142) 



2a a 

H At, iirr - 7^ — 7^ 



Q ■ r 
drgrrdrQ — 87ra{dr(pY 



4a , " n 

d^fi-T H Orgrr 

rgx rgrr 



dtKr 



^ ^ a o9 2/3 3a ^ a 
- T^d^.gr + —Kt d^gr + - 

Z,grr ^grr ^ 



gr 

grr 



a 



+ KrrKT 

grr 



a 



-{drgrf 



-drgrdrQ '^grdrQ, 



dtf 



'^grrgr'^ '^grrQ rgrrQ^ 

A3 4 2a 4a 4a 

I3drf + fdrl3 + 4 - -drl3 + i^T^.^.. + — ri^T^.Q Kt 

r grgrr gxQ rgr 

2a Aa 2a 
H K-rr, Orgx H Krr H — TrOrgrKT — IGavrllar.^, 

gxgrr rgrr gT 

I3dr4> — all, 



pdrii dl 

grr 



2a 



■drgrdr 



2a 

rg^r 



a 



Qgr 



-drQdrcl) + airn, 



grrgT 

where K = K„/ grr + 2Kt/ gr- 

The Hamiltonian and momentum constraints and the constraint defining / are 

1 1 ^ drgrr ^drgT , , 2KrrKT 



(143) 

(144) 

(145) 
(146) 



C ^ - ^^^^ + ^^^^^^^^^ 



grrgr '^grr9T T^gx T^g 

^grrg^ V 9rr 

9t \ gl gxgr 



rgrr 



rgrrgT gr 



+ 



grrgr 



(147) 



0, 



^ , 2 fKrr Kt 

drgr + 

r \grr gr 



SnUdr 



gr 



(148) 
(149) 



System (jl4Up - (jl46|) is symmetric hyperbohc with a symmetric hyperbohc constraint 
evolution system. The characteristic speeds and variables are 



f3± 



(3± 



a 



/ grr 

a 

I grr 

a 



/gr 



W 



Wn 



wz 



drgr =F 2y^i^T, 

- drgrr ~F 2-^ grr Krr '^grrfi 

■ drcj) =F V^n. 



The characteristic constraint variables are 

c' = g, 

C Cr -p "sfgrrC 



2Kt drQT 



gT 



grrgr r^gr 



(150) 
(151) 

(152) 

(153) 

(154) 
(155) 
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Note that in spherical symmetry, some of the highest derivative terms that are present 
generically in 3D cancel, both in the main evolution equations and in the implied 
constraint evolution. In particular, the right-hand side of C and Ci does not contain 
Qijk terms. For this reason Q appears undifferentiated in C^, whereas in 3D its first 
derivatives would appear. 

If we introduce an artificial outer boundary and assume that < (3 < oij ^Jg^r at the 
outer boundary, we have four incoming characteristic variables for the main system and 
two incoming characteristic constraints. The constraint preserving boundary conditions 
are 

Wtr = gtr^ (156) 

< = 9t. (157) 
6^ = 0, (158) 

Cr - y/g^C = 0. (159) 

where g^^, and g^ are two freely specifiable functions compatible with the initial data. 

We generalize the discrete boundary conditions introduced in the first part of the 
paper in the following way. We extrapolate all fields using fifth order extrapolation and 
populate grr, and at the two ghost points, r^+i and r7v+2, by solving 

Do (l - ^D+D^ grr - 2^.Krr - 2(?,,/ = (160) 



6 

h''D^grr\N+2 = 0, (161) 

Do (l - ^D+D_] - = g^, (162) 



6 

/^'/^' = 0, (163) 
where (jl60|) and (jl62p are evaluated at r^v- We then solve (jl59|) for d'^gx and use its 
fourth order accurate approximation, combined with sixth order extrapolation at riv+2, 
to populate gT at the ghost zones (we are able to compute the first derivative of gr 
from the extrapolation and the first derivatives of grr and (j) from the two maximally 
dissipative boundary conditions above). Finally, we impose ^ = by giving data to / 
at r^v using 

2 / \ 

f = Do 1 - ^D+D_ gr. (164) 

We evolve Schwarzschild space-time in standard Kerr-Schild coordinates, obtaining 
the initial data and the fixed lapse and shift from the exact solution. We introduce a 
space-like (outflow) boundary to excise the singularity. At this boundary we use fifth 
order extrapolation for the extrinsic curvature, / and 11 and sixth order for the 3-metric 
components and 0. In order to obtain a non-vacuum solution, we inject a scalar field 
pulse through the outer boundary by using 

= Asin«(t-to), to<t<to + 7r (165) 
where A = 0.02 and to = 0.1, and monitor the L2-norm of the Hamiltonian constraint at 
times t = 8M, after the pulse has completely entered the domain, and t = 16M, after 
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the pulse has fallen into the black hole. See table |21 We obtain similar convergence 
rates for Cr and Q. 

The domain extends from r = 1.8M and r = 11. 8M. We set M = 1, use a 
Courant factor of 0.1 and a dissipation parameter cr = 0.05. Dissipation is switched off 
at the last few grid points near the outer boundary, as we have observed a numerical 
instability when Kreiss-Oliger dissipation is used in combination with fifth and sixth 
order extrapolation. On the other hand, according to our tests a lower order of 
extrapolation only gives third order convergence for the contraints. 

Table 2. We give initial data corresponding to a Schwarzschild black hole in Kerr- 
Schild coordinates, inject a scalar field pulse through the outer boundary and monitor 
the constraints. Here we give the discrete L2-norm and convergence rate for the 
Hamiltonian constraint. 



Errors and convergence rates 



N 


h 




Q 


(a) t = 


8M 






100 


3.6882410" 


-4 


3.6358 


200 


2.9670410" 


-5 


4.4761 


400 


1.33308 10" 


-6 


4.3550 


800 


6.51396 10" 


-8 




(b)t = 


16M 






100 


9.88525 10" 


-5 


3.9867 


200 


6.2352710" 


-6 


4.3909 


400 


2.97205 10" 


-7 


4.4295 


800 


1.3792410" 


-8 





7. Conclusions 

We have constructed second and fourth order accurate approximations of maximally 
dissipative boundary conditions for the shifted wave equation and have shown strong 
stability of the semi-discrete scheme using the Laplace transform method. Two 
important steps were obtaining an appropriate discrete reduction to first order, and 
verifying that the solution of the problem with trivial initial data and no forcing term 
can be estimated in terms of the boundary data (the Kreiss condition). The Kreiss 
condition was verified by plotting the function N{s) to show that it is bounded. We 
have also shown numerically that the fully discrete schemes obtained by integrating our 
semi-discrete systems with a fourth-order Runge-Kutta time integrator are second or 
fourth-order accurate. 

Whereas semi-discrete approximations of one-dimensional first-order hyperbolic 
systems with constant coefficients can be transformed by a change of variables into a 
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set of advection equations, one for each characteristic variable, the same is not true for 
second-order systems§. Therefore numerical methods for second-order in space systems 
need to be developed independently. Our results for the second-order in space wave 
equation suggest a general heuristic prescription for discretising boundary conditions in 
general second-order in space hyperbolic systems: Populate all necessary ghost zones 
by extrapolation, combined with centered difference approximations of the continuum 
boundary conditions at the boundary point. We have given a simple example, the NOR 
system with scalar field matter in spherical symmetry, where this prescription gives a 
stable fourth-order accurate scheme. 

The interior schemes discussed in this paper always use a minimal width centered 
discretization, even in the outflow case {jS > 1). Despite reports on the benefits of 
the upwind (one-sided) discretization of the shift terms [211122], we find that, at least 
in the linear constant coefficient case, this is not necessary. Note that in contrast to 
an upwind scheme our semi-discrete central scheme is non-dissipative, as the Cauchy 
problem (without boundaries) admits a conserved energy. 
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Appendix A. The discrete energy method and summation by parts 

To derive continuum energy estimates one uses integration by parts to generate 
boundary terms which can be controlled by imposing suitable boundary conditions. To 
obtain similar estimate for the semi-discrete problem one can use difference operators 
satisfying the summation by parts rule jTH]. Discretisations of d^/dx^ with that property 
have been constructed in [TK] . 

In this Appendix, we attempt to use summation by parts methods to construct 
stable (and sufficiently accurate) schemes for the shifted wave equation. We are not 
successful, and try to explain why. 

We adopt the notation of reference ^H] and definitions used there, and write grid 
functions as column vectors, so that for two grid functions Uj and Vj 

N 

{u,v) = ''^^UjVj = u^v. (A.l) 

In contrast to the body of this paper, we assume that there are two boundaries in x, 
but the case with only one boundary can be obtained trivially by setting N ^ oo 

§ As an example consider the system ^(pj/At — lij, dllj/dt = D+D^ipj, whose characteristic variables 
in Fourier space are II± 2i//i sin(^/2)(/). It is not possible to translate these variables back into physical 
space. 
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and ignoring all terms arising at the right boundary in what follows. A more general 
inner product on grid functions can be characterised by a positive definite symmetric 
(A^ + 1) X (A^ + 1) matrix H, 

(u, v)h = u^Hv, = H>0. (A.2) 

Its continuum limit should be the L2 inner product. In a similar way, discrete derivative 
operators can be written as matrices. 

In this notation, the semi-discrete shifted wave equation can be written as 

^0 =/3A0 + n, 

= (3DiU + D2(f), 

at 

where Di and Di are approximations to d/dx and D2 is an approximation to d^/dx"^. 
Similarly, we can write the discrete energy as 

E = U^HU + (j)^A(t). (A.3) 

where = A and (j)^A(f) represents J {dx4>)'^dx. In particular, it should have the 
positivity properties 

^^A(f) > for all (p, (A.4) 
4>^A(f) = if and only if D+^j = (A.5) 

We assume j3 > 1 and look for (sufficiently accurate) difference operators and 
matrices H and A that give a discrete energy estimate. Taking a time derivative of the 
discrete energy ()A.3|) . we have 

—E = 2pm^HDin + (p'^ADicp) (A.6) 
dt 

+ 2U^HD2(f) + 2U^A(f). 

The key point is that the requirement that this expression be a pure boundary term leads 
to only one summation by parts condition for (3 = 0, but to three separate conditions 
for /3 ^ 0, namely 

2U^HDiU = U^\^, (A.7) 
U''HD2<P = -U''A<j)+{US<j))\^, (A.8) 
2(j)^ADi(f) = (A.9) 

where S(j) is an approximation of dx4> at the boundaries. 

We now assume that (in the terminology of ^Hl) Di is a first derivative summation 
by parts operator, meaning that 

HDi + {HDif = B, (A.IO) 

where B = diag(— 1, 0, . . . , 0, 1). We also assume that D2 is a symmetric second 
derivative summation by parts operator, i.e., it satisfies 

HD2={-A + BS). (A.ll) 
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These two assumptions guarantee that ()A.7|) and ()A.8|) are satisfied. The third 
condition, equation ()A.9|) . gives 

ADi + {ADif = S^BS + M, (A.12) 

where M = <0. 

Assuming a minimal width second order accurate approximation in the interior, we 
have that ^H] 

H = /idiag(l/2,l,..., 1,1/2), (A.13) 

A ={D+,Do,...,Do,D_), (A.14) 

D2 =iDl,D+D.,...,D+D^,Dl), (A.15) 

BS = {D+ - h/2 D^, 0, . . . , 0, D_ - h/2 D'i), (A.16) 

and (j)^A(f) = Yl!j=o{^+^jY^- ^^^^ the matrix A satisfies the desired properties. 

We need to construct an operator Di approximating which coincides with Dq at the 
interior. We consider 

(Di0)o = + (1 - (A.17) 

(^10)1 = 6^+00 + (1 - h)D+<P^ (A.18) 

and find that if 2a = 86 = 3, condition ()A.9|) is satisfied. However, this approximation 
is not accurate enough: it is only first order convergent. 
Alternatively, one could consider the modification 

(A0)o = ^/^+0o - ^^+01, (A.19) 

(Di0)i = Do0i + ah^Dl(j)o. (A.20) 

We find that M in equation ()A.12|) is indefinite (the product of two of its eigenvalues is 
negative) unless a = 0, in which case it is semi-definite positive. 

Clearly there is an infinity of other choices to make, and it is difficult to exhaust 
all possibilities. We have made a number of attempts aimed at obtaining summation by 
parts for the shifted wave equation, but we have not been able to construct operators 
and scalar products which give second order convergent schemes. We also note that a 
direct use of the operators of Appendix C.l of [T3], i.e., using and to approximate 
first and second derivatives at the boundary, gives rise to a first order refiection from 
the boundary. With the same settings used in Section we carried out convergence 
tests for /3 = 2 to confirm this. See table lAll 



Appendix B. The advection equation 

In this Appendix we list a number of known results regarding the discrete boundary 
treatment for the advection equation and give a simple prescription for the infiow fourth 
order accurate case. 
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Table Al. Using and to approximate first and second derivatives at the 
outflow boundary gives only flrst order accuracy. Here we used (3 — 2. 

Errors and convergence rates 
T2 ~q 



50 6.5465510-1 1.0630 

100 3.1334010-1 1.0252 

200 1.5395010-1 1.0065 

400 7.66263 10-^ 



Appendix B.l. Second order accuracy 

The semi-discrete initial-boundary value problem with an outflow boundary, 

jV,=aDoVj J = 0,1,2,..., (B.l) 

hmiv^i = (B.2) 

with a > is strongly stable and second order convergent for q > 2. For an inflow 
boundary, a < 0, the scheme 

= aDoVj j = 1,2,... (B.3) 

vo = 9 (B.4) 

is strongly stable and second order convergent. 

However, this last scheme is actually strongly stable for any a, even if we impose 
data on an outflow boundary. Moreover, if the redundant boundary data is second 
order accurate with respect to the continuum solution, the scheme is also second order 
convergent. 

Appendix B.2. Fourth order accuracy 

As shown in jT^j the semi-discrete initial-boundary value problem 

f^v, = aDo (^1 - y ^+^-) J = 0, 1, 2, ... , (B.5) 

hWlv^i = hWlv^2 = (B.6) 

with g > 4 and a > is strongly stable and fourth order convergent. 
In the inflow case (a < 0) the boundary conditions 

vo = g, (B.7) 
hmlv^i = (B.8) 

lead to strong stability for g = 4 or 5 and fourth order convergence is obtained provided 
that g is fourth order accurate. Interestingly, for g > 6 the scheme is unstable. One 
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can show that for a = —1, q = 6, g = 0, the scheme admits solutions of the form e**/^, 
where sh ^ 0.093 ± 1.308i. 

Note that a set of boundary conditions which lead to strong stability for any value 
of a (including the case where boundary conditions are imposed at an outflow boundary) 
is 

v-i = g^i, vo = go (B.9) 

If, in addition, g_i = u{t, —h) + Oiji^) and g^ = u{t,0) + ©(/i"^) the scheme is fourth 
order convergent [T^ . 

Appendix C. Direct second-order treatment of the shifted wave equation 

Eliminating 11 from (0121), we obtain 

9^0 = 2pdtd,(P + (1 - P^)dl^. (C.l) 

In ^2], this equation is investigated as a toy model for the Einstein equations in the 
presence of a shift. The minimal width second-order accurate discretisation in space is 

d^<j) = 2(3dtDo(j) + (1 - (3^)D+D^<j), (C.2) 

but this is stable only for \j3\ < 1. Therefore in regions where \(3\ > 1, the discretisation 

9,20 = 2/3dtDo(j) + {D+D^ - l3^Dl)(P, (C.3) 

is used instead, with a blending between the two algorithms in a transition region. The 
authors can then construct second-order accurate and stable boundary treatments for 
both outflow and time-like boundaries. They focus on Dirichlet and Neumann boundary 
conditions. 

The system could be reduced to first order in time, second order in space form by 
introducing dtcf) as an auxiliary variable. On the level of a semi-discrete (continuous in 
time) system this change of variable is trivial. In particular, the difficulty for \I3\ > 1 and 
its resolution would be the same. The crucial difference is not that ()C.1|) is in second- 
order in time form but that 11 has been replaced by dt4> in the equivalent first-order in 
time, second-order in space system. 
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